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Abstract 

In this paper we study a system of coupled nonlinear Schrodinger equations mod- 
elling a quantum degenerate mixture of bosons and fermions. We analyze the sta- 
bility of plane waves, give precise conditions for the existence of solitons and write 
explicit solutions in the form of periodic waves. We also check that the solitons ob- 
served previously in numerical simulations of the model correspond exactly to our 
explicit solutions and see how plane waves destabilize to form periodic waves. 
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1 Introduction 



Bose-Einstein condensates (BEC) made of ultracold atomic alkali gases have 
proven to be a fertile field in the last years for the study of nonlinear matter 
waves. In this context, many types of nonlinear structures have been predicted 
to exist and/or experimentally observed. Between them we can cite dark [1,2,3] 
and bright [4,5,6] solitons, gap solitons [7], domain wall solitons [8,9], dark- 
bright solitons [10,11], stabilized solitons [12,13,14], shock waves [15,16,17], 
etc. 



One of the reasons of such fruitful interplay between theory and experiment is 
the fact that below the transition temperature for Bose-Einstein condensation 
in which all bosons occupy the same quantum state, these atomic quantum 
gases are very accurately described by the so-called Gross-Pitaevskii equation 
for the wavefunction describing the colectivity of atoms thus leading to a 
complex quantum system with a highly predictable dynamics. 
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Recently, there has been a strong interest on quantum degenerate mixtures of 
Bosons and Fermions. There are different theoretical descriptions of these mix- 
tures depending on the different weight of physical phenomena. For instance 
when fermions are dominant a model was proposed in Ref. [18] and further 
studied in [19]. Solitons in the framework of this model have been found. 
Other models coming from Thomas-Fermi approximations can be more suit- 
able to describe stationary configurations with large numbers of bosons such 
as the one used in Ref. [20]. Finally, in this paper we will study a recently pro- 
posed model [21,22,23] which includes the dynamics of both the bosonic and 
fermionic component and is appropriate for describing more general situations, 
including those with smaller and/or comparable atom numbers described by 
the equations 
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where <j>F,j{T,T),j = l,...,iV>, with r G I 3 and r being the time in physi- 
cal units are wavefunctions describing each of the N F fermions, ^(r, t) is 
the mean field wavefunction for the bosonic component, atj,i,j = B,F are 
the scattering lengths for s-wave collisions between the atoms and /j, bf = 
mBm F /(m,B+m F ). The normalizations are J R „ \(f> F j(r,T)\ 2 dr = \ip(r,r) 
Nb- Finally, the trapping potential is given by 
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In the situations considered in this paper we will study situations with uj\ <C 
0J2 = 0U3 = uj±. In that limit, corresponding to strongly transversely confined 
systems for which system becomes quasi-one dimensional. Then we can use a 
multiscale expansion of the fully three-dimensional model equations [4] to ob- 
tain a one-dimensional version of Eqs. (1), which, in a new set of adimensional 
units given by x = r/a s and t = tu±, reads 
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with The constants QbbiQbf and g FF in Eq. (3) are given by the formulae 



2 



g BB = 2a BB /a s , a, g FF = 2a FF /a s a 2 , where a = m B /m F , 

a s = \J h/ (mfiUi). The new wavefuncions ip(x,t),(pj(x,t) describe the density 
variations of the condensate along the less confined dimension (x) and have the 
same normalization as the old ones, i. e. J R \(f>j(x, t)\ 2 dx = 1,/ R \ip(x, t)\ 2 dx = 
N B 

In all the numerical simulations to be shown in this paper we have taken 
m B = 87m p and m F = 40m p corresponding to a quantum degenerate mixture 
of 40 K and 8 Rb such as the one studied experimentally in [24,25]. Also the 
trapping frequency u± = 215 Hz coming from the same references will be 
taken as typical of relevant experimental situations. 

Eqs. (3) have been studied numerically in Ref. [21] and the formation of lo- 
calized structures containing bosons and fermions has been reported in the 
particular case in which the interspecies scattering length a BF is negative, 
which is the case of the 40 K- 87 Rb mixture [24,25]. In the simpler model con- 
sisting of two interacting bosonic species a negative interspecies scattering 
length has been proven to be able to support stable solitons [26,27] 

To our knowledge there have been no theoretical analyses of Eqs. (3) sup- 
porting the existence of localized structures. In this paper, we study Eqs. (3) 
analitically and numerically to provide some insight on the formation of non- 
linear coherent structures and use it to complement the numerical results of 
Ref. [21] concerning the formation of soliton trains. 



2 Modulational instability 



2.1 Analytical results 



Let us first consider the case with V B = V F = 0, then Eqs. (3) have trivial 
solutions with complex arbitrary constant amplitudes i/Jq and (pL of the form 



tp = tpQ exp i9 B (x, t) 
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for which we get, after simple algebraic calculations that 
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for arbitrary amplitudes ip B , 0^-. 



In what follows we will study the modulational instability of these plane- 
wave solutions. This mechanism has been recognized as the main one leading 
to the formation of localized structures in the framework of similar models 
[26,31,32,33,34,35] and even playing a key role in the transition to spatiotem- 
poral chaos in related models with non-conservative terms [36]. 

Now, we considerer small amplitude and phase perturbations, of the form 
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B = e B (t) + e0?(x,t), 

e? = e$ i (t) + eB*(x,t), 
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The resulting equations for the perturbations are linear and can be solved 
easily by transforming them to Fourier space, which is equivalent to looking 

for harmonic perturbations ip B , 9 B , <pfj and 9f a , j — 1, ,N F proportional 

to e lfcx e n ( fc )*. After solving the system, we get the dependence of the growth 
rate Q(k) of the Fourier mode of the perturbation with wavenumber k which 
is given by 

(n(kf - f F {k)) (m 2 - f B (k)) = C{k)\ (7) 

where 
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Perturbations remain bounded when Re[fi(/c)] < 0. The so-called modu- 
lational instability occurs when il(k) 2 > for any k. The most unstable 
wavenumbers are those with small k and a simple analysis of Eq. (7) leads 
to the condition 

9 2 bf > 0. (9) 
Then, for any g BF ^ 0, plane waves are modulationally unstable. The physical 
implications of this condition are that plane waves tend to destabilize and 
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Fig. 1. (a) Formation of trains of localized waves by modulational instability (MI) 
of plane wave initial data perturbed with a small noise of amplitude 0.015. Param- 
eter values are gss = 0.0056 and gBF = —0.0514, (corresponding to a physical 
situation witha^s = 98ao,aBF = — 410ao, with ao the Bohr's radius) Nb = 300 
and Np = 100, respectively . The solid green line corresponds to the Bosonic com- 
ponent, whereas the dashed red line identifies the Fermionic component, (b) 3D 
surface plot showing the temporal development of the MI into a train of localized 
wavepackets. 

form domains either with overlapping (gsF < 0) or separated components 
{qbf > 0). This is the reason why, in the numerical simulations of Ref. [21] 
the broad background component splits into a series of domains which we 
will identify in the following subsection with soliton solutions of the model 
equations. 

The fact that this model does not have modulationally stable solutions is 
interesting and different from the situation which happens in multicomponent 
Bose-Einstein condensates in which there is a typically a regime in which 
coupled plane waves are stable [9,26]. 

2.2 Numerical results 

To verify our results and get more insight on the development of the instability 
we have solved numerically Eqs. (3) with slightly perturbed plane wave initial 
data with periodic boundary conditions on the spatial domain x G [—30, 30] 
for different parameter values. We have observe this instability in all of our 
simulations, which confirms our analytical analysis. An example is shown in 
Fig. 1 where it is also manifest a clear tendency to form trains of localized 
wavepackets (or periodic waves) 

To observe the generation of a few solitons we have taken as initial data the 
ground state confined in a parabolic trap of unit frequency (in our adimen- 
sional units) and weak intercomponent interaction. In this initial situation we 
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Fig. 2. (a) Density profile of the ground state of a harmonically confined 
quasi one-dimensional Bose-Fermi mixture for Nb = 200 bosons and Np = 99 
fermions and gsF = —0.0038 (corresponding to an interspecies scattering length 
of clbf = — 30ao). (b) Formation of three vector solitons for gsF = —0.1253 
(corresponding to cibf = — 1000ao). (c) Formation of a single vector soliton for 
gBF = —0.0627 (for clef = — 500ao)- The other parameters are taken as in Fig. 1 



switch off the trapping potential and simultaneously increase the attractive 
interaction by changing the inter-species scattering length to a value favoring 
the formation of solitons. 

In Fig. 2 we plot the density profile of the ground state [Fig. 2 (a)] and 
the outcome after a reduction of the interspecies scattering length to clbf = 
— 1000ao, leading to three localized structures (solitons), [Fig. 2. (b)] and to 
dBF = — 500ao leading to the formation of a single soliton. [see Fig. 2 (c)]. 



3 Vector soliton solutions 



It is easy to write single-soliton solutions to Eqs. (3). We present a few ones 
here in the static case, but they can be easily modified to obtain travelling 
ones. 



First, solutions of the form of coupled bright solitons exist when g BF < 0, and 
< 9BB m B + m F \9bf\, and their explicit forms are 



ip = Aexp (lust) — -, , (10a) 

cosh rjx 

(p j = A j exp(iu F t) — - — , j = l,...,N F . (10b) 
cosh 7]X 

with parameters given by 
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This regime corresponds to the numerical experiments of Ref. [21]. Our equa- 
tions (11) provide analytical profiles for the shapes of the solitons found in 
that paper. We have checked with our numerical solutions (e.g. those pre- 
sented in Fig. 2) that the outcome of the modulational instability process 
corresponds allways to solitons very well described by our equations. Thus, 
although one could think of other more general vector soliton solutions such 
as those discussed in Ref. [37] for two-component systems they do not seem 
to arise naturally by the MI mechanism in our model. 

Although the previous case is probably one of the most interesting situations 
from the point of view of present experimental capabilities we also write here 
explicitly the form of vector solitons solutions for other scenarios for complete- 
ness. For instance bright (bosons)-dark (fermions) solutions to Eqs. (3) can be 
constructed when g BB < 0, g BF < 0, m F \g BF \ < m B \g BB \, and are given by 



ip = Aexp (iu B t) — , (12a) 

cosh T]X 

4>j = Aj exp {iujjt + i£jx) [£j + irj tsuahr/x] , j — 1, N F , (12b) 
where 
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uj b ■■ 
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Other interesting types of vector solitons of Eqs. (3) are coupled dark (bosons)- 
bright (fermions) solutions, given explictly by 
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if) = Aexp {lust + i^bx) [£b + irjtanhrjx] , (14a) 

<f)j = Aj exp (iu F t) — \ — , j = 1, N F . (14b) 
cosh T)X 

This solution exists when qbb > 0, > 0, and m F gB F < rnsgsB and the 
parameters of the solution are given by 
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Finally, explicit coupled dark soliton solutions having the form 



ip = A exp (iu B t + j^i) [£b + if] tanhr/x] , (16a) 
4>j = Aj exp {iujjt + i^jx) + i?7 tanhryx] j = 1, N F , (16b) 

can also be constructed when < gs F , msgBB < mFgBF, the parameters of 
the solution being 
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This shows that many different types of vector soliton solutions can be ob- 
tained in the framework of this model. 
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4 Quasiperiodic solutions. 



4-1 General framework 



The multipeaked structures arising in modulational instability of plane waves 
can be seen as a soliton train but are most likely to be related to the so 
called periodic solutions in the form of Jacobi elliptic functions. This class of 
solutions arise in many different related nonlinear wave equations [28,29,30] 
and we will write here some of them for our model equations (3). 

We will look for solutions of Eqs. (3) of the form 



^ = ■01 (18a) 
^■ = A^ 2 , j = l,...,N F , (18b) 



then, we get 



idtipj = --— d xx ijj + I g jk |^fc| 2 I ipj, j = 1,2. (19) 
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with mi = m B ,m 2 = m F and 
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This means that we are constructing solutions for the boson-fermion coupled 
system from those of a simpler two-component case given by Eqs. (19) and 
(20). In particular we will look for solutions of the form 

ipj = Cj exp (iujt) , fj (£x) (21) 

where fj is one of the cnoidal functions fj = sn, cn, dn satisfying the ODE 

/r (", • '•,./;) ./:,• (22) 

and any two different of them are related by 

Pifi+P2fl = l, (23) 
where Uj, Vj and pj are real numbers. 
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4-2 Cases with fx = f 2 . 



First we consider the situation in which ipi and ip 2 ar e proportional, i.e. f\ = 
f 2 = f — sn, cn, dn. Then, after substituting (21) into Eq. (19) one gets that 
the dispersion relation reads 



(cf) 


C2 ( 




2gi292i y 



"* - (24) 

and the amplitudes are given by the equations 

Q\ 2 lm 2 i 

' (25) 

where u = u\ = u 2 and v — v\ — v 2 . There are different subcases depending 
on the particular choice of /. Taking / = sn(£a;, k) we get 

u = -(l + k 2 ),v = 2k 2 . (26a) 

This solution exists if 

Qbf > 0, m B g BB < m F g BF . (26b) 

Next we consider / = cn(£x, k) for which we get 

u = 2k 2 — l,v — —2k 2 , (27a) 

which exists provided 

Qbf < 0, m B g BB > m F g BF . (27b) 

The last possibility with fi = f 2 corresponds to / = dn(£x, k), its dispersion 
relation being ruled by 

u = 2 — k 2 , v = —2. (28a) 

This solution exists if 

9bf < 0, m B g BB > m F g BF . (28b) 
4-3 Cases with f\ ^ f 2 . 



In the more general situation j\ ^ f 2 we have a larger variety of combinations. 
In what follows we write the explicit formulae for those situations. First, after 
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substituting (21) into Eq. (19) we get the dispersion relations, 



e 



2m i 



Uj + 



9jjC 



while the amplitudes are given by 

9uPi/P2\ (vi/mi 
\92iP2/p\ 9u J \V2/m 2/ 
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(29a) 



(29b) 



Now our task is to analyze the different possible combinations of sn-, cn- and 
dn-functions. 



4-3.1 Case f\ = sn(^x, k), f 2 = cn(£x, k). 

The parameters in the dispersion relation are given by 

vi = 2k 2 , v 2 = -2k 2 , pi= p 2 = 1. 
The amplitudes of the solutions satisfy 

/ „ /™ \ 



CI 



2U2 



912921 



9i2 /m 2 
\9ii/iri 2 - g 21 /m 1 J 



This solution exists if 



g B F > 0, m B g BB > m F g BF . 



(30a) 



(30b) 



(30c) 



4-3.2 Case f\ = sn((,x, k), f 2 = dn^x, k). 

The parameters in the dispersion relation are given by 

vi = 2k 2 , v 2 = -2, pi = k 2 , p 2 
The amplitudes of the solutions satisfy 

( 



1. 



CI 



912921 



gi 2 k 2 /m 2 
\gu/m 2 - g 2 \lm X/ 



This solution exists if 



(31a) 



(31b) 



g BF > 0, m B g BB > m F g BF . 



(31c) 
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4-3.3 Case f\ = cn(£x, k), f 2 = sn^x, k). 

The parameters in the dispersion relation are given by 

v 1 = -2k 2 , v 2 = 2k 2 , Pl = p 2 = 1. 
The amplitudes of the solutions satisfy 

\-gu/m 2 + gixjmi) 
This solution exists if 

9bf<0, 9bb<0, m B \g BB \ > m F \g BF \ 

4-3.4 Case f\ = cn{^x,k), f 2 = dn{^x,k). 

The parameters in the dispersion relation are given by 

k 2 1 



6) 


_ ek 2 ( 







Vl = -2k\ v 2 = -2, 



Pi 



1 - k 2 ' 



P2 



1 - 



The amplitudes of the solutions satisfy 
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c 2 i fj^fj-^ 



;2 ( -g 12 k 2 /m 2 ^ 
K gu/m 2 - W m iy 



This solution exists if 



9bf < 0, m B g BB > -m F \g BF \ 



4-3.5 Case f\ = dn(£x, k), f 2 = sn(^x, k) . 

The parameters in the dispersion relation are given by 

v 1 = -2, v 2 = 2k 2 , p 1 = 1, p 2 = k 2 . 
The amplitudes of the solutions satisfy 



( 



-912/ m 2 



\ 



K k 2 {g 2 i/m 1 - g n /m 2 ) J 



912921 

/ 

This solution exists if 

9bf < 0, gBB<0, m B \g BB \> m F \g BF \ 
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m B 9BB > Ttip \qbf\ 


cn 
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Table 1 

Summary of exact solutions in form of coupled cnoidal waves to Eqs. (3) with (18). 



4-3.6 Case f\ = dn(^x,k), = cn{^x,k). 

The parameters in the dispersion relation are given by 

1 k 2 



ui = -2, v 2 = -2k', 



Pi 



1-fc 2 ' 



P2 



1-k 2 



The amplitudes of the solutions satisfy 
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- e { 
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This solution exists if 



9bf < 0, m B g BB > -m F \g BF \ . 



(35a) 



(35b) 



(35c) 



4-4 Summary 



Table 1 lists the different possible combinations found in the previous subsec- 
tions 
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5 Conclusions 



In this paper we have studied a mean field model proposed to model the 
dynamics of Boson-Fermion mixtures. We have found that plane waves are 
allways modulationally unstable in this model and have shown by means of 
numerical simulations how their destabilization gives rise to structures sim- 
ilar to periodic waves. We have also studied numerically the destabilization 
of the ground state of the system and the generation of individual solitons. 
The explicit form of these solitons have been obtained for the case at hand 
and for other possible combinations. Many of these soliton solutions of our 
model equations can be observed in future experiments with Boson-Fermion 
mixtures. Finally, many explicit solutions in the form of periodic waves (sn,cn 
and dn functions) have been constructed. 
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